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1.  INTRODUCTION 


The  basic  problem  of  periodic  interpolation  can  be  explained  as 
follows : 

-  given  a  finite  set  of  distinct  points  t ^  (  ...tt  of  the 

interval  and  a  set  of  real  scalars  / 

-  construct  a  iff  -  periodic  function  s  belonging  to  the  linear 
variety  V  .  of  all  interpolants  of  a  suitable  space 

1/  .  /  /  £  #  I  fuj  .  A/}.  (1.1) 


For  practically  relevant  candidates  of  solution  additional  information  is 
clearly  needed,  in  particular  polynomial  precision  and  smoothness.  This 
essentially  amounts  to  restricting  the  set  of  ?€  -  interpolants  by 

constraints  involving  a  suitable  “energy"  norm. 

In  this  connection  a  quadratic  (semi-)  norm  generated  by  the  integral 

{/  I  (1.2 


has  proved  most  efficient,  where  the  differential  operator 


f) ...  ({ 


afx/ 


(1.3) 


annihilates  all  trigonometric  polynomials 

m 

S(x)  *  (  a  k  CO$(kx)  be  sin  (k*)) 

it  »  O 


2 


of  order  m. 

A  proper  setting  for  the  problem  of  minimizing  the  (quadratic,) 
functional  (1.2)  under  interpolatory  constraints  is  provided  by  the 
class  7C  =  irj  Of  zr  -  periodic  functions,  whose 

(distributional)  derivatives  up  to  the  (2m  +  2)  -th  order  are  square- 
integrable  on  the  interval  £o,  Jtt]  . 

The  solution  s  of  the  minimization  problem 


J.T 

f  13^  sM^dx 


O 


cn/  J  13^  ^(x)^  dx 

4***  ° 


(1.5) 


is  called  (optimal)  trigonometric  spline  interpolant  in  X  (  ”*  ^[OtZir] . 

The  trigonometric  spline  interpolant  is  essentially  given  by  the 
following  properties  (cf.  Schoenberg  (1964)): 

(i)  s  is  ifr -  periodic  and  continuous 

(ii)  s  is  infinitely  differentiable  for  all  points  -t  e  Co,2v) 

+  tk  f  k  ,N  with  SCO  =  O, 

i.e.:  s  reduces  to  a  trigonometric  polynomial  of  order  m  for 
all  points  •£  €  [O,  Sir)  with  ****,*-*,...,  v. 

( i i  i )  s  (**)  *  yk  for  k  *  -f,  . .  .  ,  V  . 


Spline  interpolation  turns  out  to  be  a  most  adaptable  method  to 
data  for  (global)  interpolation,  and  the  (semi-)  norm  (1.2)  is  a  natural 
setting  to  maintain  tne  flexibility  of  piecewise  trigonometric  polynomials 
while  at  the  same  time  achieving  some  degree  of  global  smoothness.  Of 
particular  usefulness  is  the  spline  interpolant  correspondi ng  to  the 
Quadratic  integral  mean  of  the  (linearized)  curvature  (et/tfx)*' f 
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/ 


/3„ 


(1.6) 


since  the  quantity  (1.6)  may  be  physically  interpreted  as  the  potential 
energy  of  alstatically  deflected) thin  beam  (which  indeed  is  proportional 
to  the  integral  of  the  square  of  the  (1 inearized)  curvature  of  the 
elastica  of  the  beam)  (cf.  Moritz  (1978)). 

Trigonometric  splines  may  be  interpreted  as  the  spline  functions 
for  the  (unit)  circle,  i.e.  trigonometric  splines  are  the  two-dimensional 
analogues  of  the  spherical  splines  discussed  in  this  paper. 

o 

Roughly  speaking,  the  spherical  i nterpol ation  problem  to  be  of 
importance  for  geodetically  relevant  purposes  can  be  formulated  as  follows: 

-  given  a  finite  set  of  distinct  points  ( stations )  f  •  •  • ,  y ^  of 
the  (unit)  sphere  £2  and  a  set  X,  ,  *  •  •  ,  Y,j  of  real  scalars 
(observations  or  measurements) 

-  find  a  (smooth)  function  5  £2  — —  ^  belonging  to 

the  linear  variety  of  all  interpolants  of  a  suitable  space  <a? 


- 


=  {  4  e  %  i  l  ( )  «  yk  ,  **  *  'f, 


A 1 


(1.7) 


In  order  to  achieve  uniqueness  it  likewise  seems  quite  natural  to  look  for 
an  interpolant  minimizing  an  appropriate  quadratic  norm  in  such  a  way  that 
additional  assumptions  concerning  polynomial  precision  and  smoothness 
again  are  satisfied.  This  can  be  done  in  the  same  way  as  in  the  trigono¬ 
metric  case  by  using  a  differential  operator  annihilating  all  polynomials, 
i.e.  spherical  harmonics  of  order  m  or  less.  Observing  the  fact  that  the 


* 
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spherical  harmonics  Sn  of  order  n  are  the  everywhere  on  the  (unit) 
sphere  regular  eigenfunctions  of  the  Beltrami  operator  A*  corresponding 
to  the  eigenvalues  A  ^  -  «  (  *n.  +  >i  )  (  i  .e. 

(A-  *  XJ  (f)  =  O  (1.8) 


on  the  unit  sphere  Q  for  all  spherical  harmonics  5  of  order  n,  the 
(simplest)  operator  playing  the  same  role  as  the  operator  (1.3)  is  given  by 
the  product 

(A')  -  +  ...  (  A*  +  Am) .  (1.9) 

rry 


A  proper  setting  for  optimal  spherical  interpolation  is  provided  by  the 
(Sobolev)  space  —  7C  (Xn*  *X)(  )  of  functions  whose  (Beltrami) 

derivatives  up  to  (A*)  4  (in  the  distributional  sense)  are  square- 

integrable  on  the  sphere  £1  (cf.  Chapt.  8).  Equipped  with  the  (semi-) 
norm  f  •  generated  by 


the  linear  space  ^  (  SZ)  is  a  (semi-)  Hilbert  function 

space  of  continuous  functions  on  the  (unit)  sphere  Q>  .  Therefore,  in 
comparison  with  (1.5),  the(optimal)  spherical  interpolation  problem  to 
be  analyzed  in  this  paper  can  be  formulated  as  follows: 

-  find  a  function  s  e  such  that 

Ii(A*.)  -inf  i  I  (Al\)  . 

a  '■  fe  vu  a 


(1.11) 


5 

As  described  in  detail  the  variational  problem  (1.11)  is  well-posed  in 
the  sense  that  its  solution  exists,  is  unique,  and  depends  continuously  on 
the  data  Xt  ,  •  .  •  ,  >,v  •  Basic  tool  is  the  theory  of  Green's  functions 
of  the  sphere  with  respect  to  the  operator  (1.9). 

The  spherical  spline  interpolant  is  essentially  given  by  the  following 
properties: 

(i)  s  is  continuous  on  £2 

(ii)  s  is  infinitely  differentiable  for  all  points  ^  e  S2  ,  7$, 

it  -  '1  with  (A’U,)4  .  .  .  (  +  =  0,  i.e. 

s  reduces  to  a  polynomial  of  order  2=  m  for  all 
points  £  €  52,  with  jr  7-  7k  t  &  =  -/,•••,  /v' 

(ill)  5  ( 7h )  ~  X*  for  k  m  4 1  ;  A/  . 


Spherical  spline  functions  (s.s.f.)  have  the  following  attractive 
features: 

The  linear  space  of  all  s.s.f.  (of  order  m)  is  finite  dimensional; 
s.s.f.  are  relatively  easy  to  manipulate  and  compute;  various  matrices 
arising  in  interpolation  and  approximation  problems  are  nonsingular, 
s.s.f.  do  not  tend  to  produce  approximations  having  severe  oscillations 
and  undulations. 

Thus  s.s.f.  seem  to  be  best  suited  for  macro  -  and  micro  modeling  of 
earth’s  gravitational  field.  Moreover,  a  variety  of  problems  of  "optimal  1 
integration  on  the  sphere  leads  to  spherical  spline  ^unctions.  In 
addition,  the  whole  solution  process  can  be  made  surprisingly  simple  and 
efficient  for  numerical  computation.  Indeed,  the  computational  scheme  for 
solving  the  interpolation  problem  can  be  described  in  a  recursive  form, 
based  on  a  combination  of  generalized  Lagrange  and  Newton  formula.  This 
means  that  the  method  is  constructed  so  as  to  have  the  so  -  called 
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permanence  property:  the  transition  from  the  solution  with  respect  to 
N  data  to  the  solution  with  respect  to  (N  +  1)  -  data  necessitates  merely 
the  addition  of  one  more  term,  all  the  terms  obtained  formerly  remaining 
unchanged. 


The  price  to  be  paid  for  the  convenience  of  the  permanence  property 
is  a  bi orthonormal izati on  process. 


a 
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2.  DEFINITIONS  AND  NOTATIONS 


denotes  the  three-dimensional  real  Euclidean  space.  We  consistently 
write  x  ,  y,  i(  ...  for  the  elements  of 


Let  e.  .  e*  ,  be 

the  canonical 

/a  \ 

f° 

'  iS  ' 

II 

V  0 

>  e  i  ~ 


(2.1) 


(2.2) 


In  comDonents  we  have/for  elements  x  (  g  g 

*  =  z,  e,  +  **ez  +  **  e3  • 

0  3 

The  inner  product  of  x  (  j  e  f  is  the  number 

*  a  =  X*I  «*  £3  / 

—  3 

the  norm  of  x  £  -"O  is  the  (nonnegative)  number 

>xl  =  J  X  ■  X  .  (2.4) 

For  all  elements  A  e  ,  different  from  the  origin,  we  have  the 
representation 


(2.3) 


>  =  v 


■r  =  /  x  1  —  /  X,*"  +•  x.  +■  x  , 


)  2 . 5) 


where  is  the  uniquely  determined  directional  unit  vector  of  the 

element  x  e  /j?3. 

The  unit  sphere  in  7R  will  be  called  £ 2  .  The  total  surface  of 
Q,  wi II  be  denoted  by  : 
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&  -  J  Cttjy  *  V* »  .  (do :  surface  element). 

Q 

The  rectangular  coordinates  (2.1)  are  related  to  the  polar  coordinates  (2.5) 
by  the  equations 

x  -  *  if  ,  ^=-£e3-r  J  a-  t2-  (ax(peA  +  Sm  <pej  (2  6, 

■6  *  Coi  ■&  ,  0  T  t  O  &  <p  c  -?7T, 


.e. : 


x„  =  -r  /'f  -  f*  cos  £2? 

=  V  J  /I  -  -L*  Sin  (p 


X 


3 


T  "t 


(2.7) 


(■6s:  polar  distance,  p  :  geocentric  longitude).  In  terms  of  the  coordinates 
(2.6)  the  Laplace  -  operator 


/* 

v  3  x. 


takes  the  form 


/  2  ' 


A  H 

■*"  3-r 


(2.8) 


(2.9) 


A*  denotes  the  (Laplace  -  )  Beltrami  -  operator  of  the  unit  sphere 


^  (A 


f 


(2.10) 


5  e 


(  t 


)  ^ 


-i-t 


(2.  i 


X 


c 


9 


3.  SPHERICAL  HARMONICS 

The  (Laplace  -)  spherical  harmonicsS*  of  order  n  are  defined  as  the 
everywhere  on  the  unit  sphere  S2  infinitely  di fferentiable  eigenfunctions 
of  the  Beltrami  differential  equation 

(A"f  *  An  )  (f)  -  O  (3.1) 

corresponding  to  the  eigenvalues 

n  *  '•)  ■  (3.2) 

The  set  of  all  eigenvalues  is  the  spectrum 

=  /  Xn  =  n  (n+  -i)  |  n  -  O,  ^ I,  ■■.}  .  (3.3) 

As  is  well-known,  the  functions  H n  given  by 

=  -r*  (f)  ,  3.4) 

are  polynomials  in  rectangular  coordinates  which  satisfy  the  Laplace 
equation 


(x)  =  0  (3.5) 

and  are  homogeneous  of  degree  n.  Conversely,  every  homogeneous  harmonic 
polynomial  of  degree  n  restricted  to  the  unit  sphere  £2  is  a  spherical 
harmonic  of  order  n. 

The  Legendre  polynomials  ^  given  by 


Cn/i] 


.s  ■  o 


,  .  s  C  in-ls)\ 

(-  a) - 

2n(n-2s)'.  (n-s) ■!  s! 


t 


n  -  2s 


(3.6) 
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are  the  only  everywhere  on  the  interval  £--1,  -ij  infinitely  differentiable 
eigenfunctions  of  the  Legendre  differential  equation 


^n]  ,  -i  €  [- a,  a ]  /  (37) 


which  in  ~t  ~  ^  satisfy  ('O  —  ^ 

Apart  from  a  constant  factor,  the  Legendre  polynomials  are  the  only 
spherical  harmonics,  which  are  invariant  under  orthogonal  transformations 
with  the  "north-pole"  e3  as  fixed  point. 

o 

Spherical  harmonics  of  different  order  are  orthogonal  in  the  sense 
of  the  •£ 1  -  inner  product 


(&n  ■ 


fS«(V  fa  m  °  ■ 

SI 

( n  +■  m ) 


The  linear  space  ^  of  all  spherical  harmonics  of  order  n  has  the 
dimension 

dim  (fn  )  =  Zn  +  d  .  (3.9) 

In  other  words,  there  must  be  n  +  'f  linearly  independent  spherical 
harmonics 


-S 


in*  A 


(3.10' 
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We  assume  this  system  to  be  orthonormal i zed  in  the  sense  of  the  -inner 
product 


SI 


=  <£*,  <£h 


(3.11) 


(  ‘Cnt'-  Kronecker  symbol) 


denotes  the  linear  space  of  all  functions  in  three  variables  which 
restricted  to  the  unit  sphere  72  may  be  represented  by  homogeneous  poly¬ 
nomials  of  degree  m  or  less.  admits  the  orthogonal  decomposition 

9„  -  4  ® 

with  respect  to  the  ^ -inner  product,  i.e.  to  every  S  €•  there 

exist  spherical  harmonics  S0  (  .  .  .  ,  S  ^  with 


5  -  21  Sr 


,  <X,  V  -  0  * 

(n  *  J  ,  O 


Hence,  the  dimension  of  P  i s  equal  to 


(3.12) 


M 


1  *^1  ,  « o 


;  3 . i3) 


For  any  two  elements  ,  the  sum 

is  invariant  under  all  orthogonal  transformations  A,  i.e. 


f.  cf,  ?> 


<*■%, A?) 


(3.14) 
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for  all  orthogonal  transformations  A.  For  fixed  ^  €  £2.  >  ( if  !  y) 

is  as  function  of  %  a  spherical  harmonic  of  order  n.  ^  (it  ? )  is 
symmetric  in  If  ,  £  and  depends  only  on  the  scalar  product  of  f  and  ^ 
Thus  it  is  clear  from  the  above  that  we  have,  apart  from  a  multiplicative 
constant  «.  ^  , 


2n+ 1 


[3.15) 


In  order  to  determine  ^  we  set  if  =  £  .  T^en  we  obtain 


3n  +  4 


I  = 


CX. 


Integration  over  S2  yields 

Jn  +  A  *  V’TT  <xf 

Therefore  we  find  the  addition  theorem 


(3.16) 


+-i 


2n*  * 

~TtT 


2  <•?•?) 


(3.18; 


In  particular,  we  have 


►n 

H 

rv«  O  )•' f 


s-v  -  6.  ^ 


2 'ft). 


(3.19) 


Let  <fi  be  an  absolutely  integrable  function  on  the  interval  £- ,  ij  , 
i  .e  :  (p  6  Then,  for  every  -Sn  £  ,  Hecxe's  formula 

(cf.  Muller  (1966),  Freeaen  (1979))  gives 


<p  (%■?)  Sn  (£)  d(j({)  - 


(3.20) 
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where 


2i r  J  <p  (-t)  'PU)  d-t . 


(3.21) 


The  notation  ddy[p)  means  the  surface  element  is  applied  to  the  in¬ 
variable. 

Hecke's  formula  establishes  the  close  connection  between  the  orthogonal 
invariance  of  the  sphere  and  the  addition  theorem. 
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4.  FUNDAMENTAL  SYSTEMS 


Unfortunately,  the  system  of  spherical  harmonics 

/  ‘  *  *  1  •  I  "^n.xn-n  I  '  I  1  <  *  '  •  /  ^rvi, 

is  not  unisolvent  on  the  unit  sphere  £2  ,  i.e.  the  matrix 

v 

(. 


I  7. ) 

( 7,) 

^,3  (?«) 


\  '*'*,*■'”+*  (t,) 


S*,«  (  7m) 

(  7m) 

S,,3  (?*) 


(7m) 


(  7m) 


(4.1) 


(4.2) 


is  not  non-degenerate  for  all  choices  of  M  distinct  points  7, ,  . . . ,  7^ 
lying  on  22  .  This  is  known  from  Haar's  theorem  (cf .  Davis  (1963)  Theorem 
2.41)).  However,  it  is  easy  to  prove  that  there  exist  systems  of  points 
7„/.../7m  having  a  non-degenerate  matrix  (4.2).  For,  it  is  certainly 
possible  to  find  a  point  ^  with  -S0/,  (7^  )  =ft  o  -  We  consider  the 
determinant 


dei 


I  s«.«<7«) 

(  7i)  S,.„  ( 


(4.3) 


As  a  function  of  !•  ,  this  determinant  cannot  be  identically  0,  for  else 
So,*  and  S. ,  would  be  linearly  dependent.  Therefore  there  is  a 
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point  >r4  such  that 


de-t 


det 


\ 

+  < 

(7-) 

the  determinant 

(*.) 

'?*.) 

\ 

(7.) 

^,-t  (  7j t) 

4.,  (?) 

(?*> 

(?)/ 

(4.4) 


(4.5) 


and  by  the  same  arguments  there  exists  a  point  tf  *  for  which  (4.5) 

is  different  from  zero.  Therefore,  by  induction,  we  can  find  a  system 
of  points  such  that  the  matrix  (4.2)  is  non-degenerate. 


A  set  of  M  points  of  the  unit  sphere 

is  called  fundamental  system  of  order  m  on  &  ,  if  the  rank  of  the  (M,M) 
-  matrix 


cc 


h 


XK) 

S-M  (?.) 

(?.*) 


So,,  (?„) 

(  7m) 

^4  (?*) 

(7j 


\ 


(4.6) 


>  +  /I 


(7m)/ 


is  equal  to  M. 
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A  set  «*  ,  .  .  .  ,  7.,  of  «  (Vn+>i)i  (distinct)  points  of  the 

unit  sphere  £2  is  called  admissible  system  of  order  m  on  ^  ,  if  the 
rank  of  the  (M,N)  -  matrix 


oc 


.  .  . 

5...  (?„) 

^  ( 70 

.  •  • 

(  7/v) 

•  •  • 

C  w 

(4.7) 


is  equal  to  M. 


(•?.) 


(V^) 


NOTE:  In  the  sequel  we  assume  that  each  admissible  system  v  ,  •  •  •/ 
of  order  m  has  as  subset  the  fundamental  system  ,  ■  ■  ■  ,  ? ^  of 

order  m  (consistinq  of  the  first  M  elements  •?..  ).  This 

-  t  /f  i  ft  A"t 

is  always  achievable  by  reordering. 


Given  a  function 

s(i)  - 


S  e 


of  the  (general )  form 


ry\  in+1 

n  n 


(4.3) 


For  an  admissible  system  ^  ,  .  . 

linear  equations 


I]  ^ 

k 


In 


c f  order  m  on 


,  the 


(4.9) 

(i^o,  .  ,  m;  y  =  V,  ..  .  .inf-i) 


are  solvable  in  the  unknowns  a,  ,  .  .  . ,  .  Thus,  ^or  every  admissible 

system  ^  f  .  .  .  t  y  v  of  order  m  on  .12  and  all  solutions  ,  .  .  .  ,  cf 
the  linear  equations  (4.9)  the  element  5  €  of  the  form  (4.3! 

can  be  expressed  as  follows: 
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HI  &H+A 


S(i !  - 

*--»  "'O  /' I  'J  ' 


Using  the  addition  theorem  we  obtain 


(4.10) 


(4.11) 


But  this  means  that  the  N  functions  given  by 

•» —  j.  +  'i 


<*-Tf 


3 


(4.12) 


z:  ^ 


span  the  space  : 


O 


-  SPan-  ((.?  *5?  ^fS,lU, . „).  (4.13) 


a 
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5.  GREEN'S  FUNCTIONS 


We  next  define  the  Green  function  of  the  unit  sphere  S2  with  respect 
to  an  operator  A*  +  €  S(Q.)(cf.  Free  den  ( 1?78/19'?9) ) . 

This  function  will  be  of  great  importance  for  the  definition  and  the 
application  of  spherical  spline  functions. 

We  begin  our  considerations  by  introducing  the  definition: 

A  function  £  ( An  ■  is  called  Green's  function  of 

the  uni  t  sphere  Sc  with  respect  to  the  operator  A*  +  .A  n  and_ 
the  parameter  ^  e  S2  if  it  satisfies  the  following  properties: 

(i)  <2  (  Xn  •  is  as  function  of  *7  ,  for  fixed  £  ,  infinitely 

differentiable  for  all  *7  e  S2  with  ?  ^  f  (differentiability 

(ii)  For  all  yeQ  with  g  ?  if 

(differential  equation) . 


(iii )  For  all  y  e  SZ  , 

|  <"V  f,?)  *  f  ?) 

is  continuously  differentiable  for  all  ^  e  Q  (characteristic 
singularity) . 

(iv)  For  all  orthogonal  transformations  A 

^  (*n;  A?'  A7]  =  $  CXn  '  Y} 

(rotational  symmetry )  . 

^  J  1  n  ;  S  t  *?)  ^  (57)  l' 7)  ”0 

uniformly  with  respect  to  a  II  ^  £  Jo 

(normalization) . 
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We  first  prove  the  unioueness  of  Green's  function  ft %): 

Denote  by  £>(\n/  f  t  y)  the  difference  of  two  Green’s  functions  satisfying 
the  properties  (i)  -  (v).  Then  we  have 

(i )  '  -3)  (■  ;  f  i  •?  )  is  as  function  of  £  ,  for  fixed  %  €  Q.  , 

infinitely  differentiable  for  all  ?  €.  5>2  with  ^  ^  ^ 

( i i )  ' For  all  y  £  ^2  with  p  ^  ^ 

(A’  +  >„)  D<K,-  ?.7)  -  0  • 

( i  i  i ) '  For  all  ?  e  £2  ,  j)(An;  f,y)  is  continuously  differentiable. 

(iv) 'For  all  transformations  A 

3^;  -  3>  (**;$,?) 

(v) '  J  5)  (\n;  f,p)  ?)  doty  =  o 

uniformly  with  respect  to  all  f  €  . 

The  properties  (i)‘  -  ( i i i ) *  show  that  3)  (  A  ^ ;  y)  is  an  everywhere  on 
the  unit  sphere  £2  infinitely  differentiable  function  satisfying  the 
differential  equation  ( i  i ) '  .  Therefore  3(X„;  f,pj  must  be  a  spherical 
harmonic  of  order  n.  Because  of  the  property  ( i v )',  3)CAn/  $/?)  depends 
only  on  the  scalar  product  if-  £  .  Consequently,  we  have 

f,?)  -  rn 

From  (v)'  we  deduce  that  =  O  .  But  this  means  that  Green's  function 
^  n ■  f t  ?)  of  the  unit  sphere  with  respect  to  the  operator 

+  A  ^  and  the  parameter  is  uniquely  determined  by  the 

defining  properties  (i)  -  (v). 

Following  Hilbert's  approach  to  the  theory  of  Green's  functions  (cf. 

Hi  1  bert  ( 1912) )  we  prove  the  exi stence  of  £  ( ^  ,  p)  by  first 
giving  an  explicit  representation  of  Green's  function  to  the  operator  A  . 


An  easy  calculation  shows  that 
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—  s&tj-  (  ^  ~  f  ■  p)  —  'I  + 

satisfies  all  the  defining  properties  (i)  -  (v)  of  Green's  function  with 
respect  to  A*  .  Hence,by  virtue  of  the  uniqueness,  we  have 


$'  0,1.7)  + 


(5.1) 


In  order  to  assure  the  existence  of  ^  ^  ^ 

we  consider  the  integral  equation 


fO,,  f,  9)  ‘ 


(5.2) 


*  is.  / 


A_ 


~XT 


?« <i-?h 


which  establishes  the  close  relation  between  Green's  function  4 f,p)  and 
the  resolvent  of  the  kernel  ^  7)  •  It  is  not  difficult 

to  see  that 


f  <  '<>,•  f,  7)  ^  O) 

a  (5.3) 

-  /  [  r  *  ^  'Pjf-r^s,  <1)^(7) 

SI  H  'V 

for  all  spherical  harmonics  -S„  £  of  order  n>0.  Tne  integral  equation 
(5.2)  therefore  has  a  solution  whicn  is  uniquely  determined  by  the  conditions 


(  J  O'n ;  %  1*1)  ^ ? )  *  C 

6? 
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for  all  ^  . 

Observing  the  characteristic  singularity  of  Green's  function  we 
obtain  applying  Green's  (surface)  identity 


^  ft  ~  ^  n. 
^7 r 


/  =  ('<-JA  x  )  SnM.  (5.5) 

£ 


Thus  the  spherical  harmonics  of  order  n,  i.e.  the  eigenfunctions  of 
the  Beltrami  -  operator  A*  with  respect  to  the  eigenvalues  =  n(n  +  -i) 
are  eigenfunctions  of  Green's  (kernel)  function  (^n. ;  %  /  7 )  in  the 
sense  of  the  integral  equation  (5.5).  Furthermore,  if  £  ^7, 

with  -  'i  &  if- y  <  -/  ,  the  kernel  ^  (Xn ;  f,  ?)  allows  the  bilinear 

expansion 


5  ft  7^ 


OQ 


-  r 

/?=o 


ty-  TV 

-4* -A. 


ZT  4„(?)  5.,/?)- 

j  ~  1 


(5.6) 


The  symbol  ^L*  denotes  that  the  sum  is  to  be  extended  over  all  nonnegative 

integers  £  with  X.  A 

*<  '  n.  ^ 

u 

Using  the  addition  theorem  we  can  rewrite  the  bilinear  expansion  of 

f(A„;  ^ 1n  the  f0rTT1 


f,7>  -  Z*  /Zr  ^(??) 

0  k*o  Ak~  *n. 


(5.7) 


According  to  the  classical  Fredholm  -  Hilbert  theory  of  linear 
integral  equations  (cf.  Hilbert  (1912))  we  define  iterated  Green's  functions 
by  the  following  convolutions 

(5.8) 

=  j  y-1  (Ac  f  ■  ■  ■ ,  A  »<-■>;  £ •  ~C,  )  ( ;  *?/ 7)  o/cJ  (  7  ) 

x7 


^ ,  ■  • ) 


f(A0'...t  A  ^  .  f ,  %)  is  called  Green's  function  of  the  unit 
sphere  with  respect  to  the  operator  (A*)^  -  .  .  (A  +  Am). 

In  analogy  to  techniques  known  in  potential  theory  it  can  be  shown 
that,  for  integers  m  3*  'l  ,  $  (Aol  - ? )  is  continuous  on  the 

whole  sphere  S2  as  function  of  y  with  £  fixed  or  as  function  cf  % 
with  ^  fixed. 


On  the  other  nand,  the  bilinear  expansion  of 


is  absolutely  and  uniformly  convergent  both  in  and  £  respectively 
and  uniformly  in  $  and  £  together.  Thus  the  representation  theorem  of 
(general ized)  Fourier  theory  yields 


VX  *>■■■,  •On)  ZT  {5*10' 


Let  in  be  an  integer  with  Then  the  derivative 


(Am+A0)  ...  (  £  <^o,  J7>) 


(5.11) 


isas  function  of  ^  for  fixed  ^  ,  a  continuous  function  on  the  unit  sphere. 
For  integers  «n?1,  the  derivative 


f&m  +K)  •  •  •  (’  +  A>"-r}  ?  *  •  •  /  A*  ;  f,?} 


(5.12' 


ocssessesas  function  of  logarithmic  singularity  in  f e  .  Furthermor 

for  elements  with  g  ys  ^  and  integers  >,  -> if 
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6.  INTEGRAL  FORMULAS 


Let  if  be  a  fixed  point  of  the  unit  sphere  £2  .  If  now  -f  is  a 
function  with  continuous  second  derivatives  on  £1  ,  tnen  for  each  sufficiently 
small  jo  >  0  Green's  surface  identity  gives 

J"  { t(xo;  f,z)[Am  fa)]  -  [a*  faOJ  f, rfl  dfa\  (6.D 

(71  *  1 

■  f  1  $(k  m  -  MA  fa°>  - 

tl-w  *  1 


Herein  ds  is  the  line  element  in  7R  3  ,  while  n  is  the  unit  vector, 
normal  to  the  curve  l'f-iz  I  -  J3  on  Q,  ,  tangential  to  Q  and 
directed  exterior  to  the  set  of  all  y  €  S2.  with  . 

In  identity  (6.1)  we  first  observe  the  differential  equation  of  Green's 
function 

/  /fy  [£>  =  J  fay)d<*(y).  (6.2) 

tf-Z  izp  't-rs*P 

in  polar  coordinates  (2.6)  the  line  element  for  If’Z1  “  f  can  be 
expressed  by 


-  jfji  ^  *  ( -<  - 1*)  a  4* . 
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Therefore,  by  virtue  of  the  logarithmic  singularity  of  Green’s  function 
^  f*,  £  )  ,  we  set  in  analogy  to  the  well-known  consider¬ 

ations  of  potential  theory,  on  passing  to  the  limit  p  — O  the  theorem 

Let  if  be  a_  fixed  point  of  the  unit  sphere  £2  .  Suppose  that 

■f  is  a  twice  continuously  differentiable  function  on  £2  .  Then 

=  ^  /  P?)  dot?)  (6.3) 

SI 

Thi  s  formul  a  compares  the  functional  value  of  a_  functi  on  /  t=  C  u Y o2J 
at  <£  22  wi th  the  mean  (integral )  value  of  /  on  the  unit 

sphere  S2 .  . 


3y  use  of  Green's  function  f  ^o,  ;  f  ,  £  )  with 

resoect  to  the  operator  (A*+A Q)(A*+A1)  we  are  able  to  generalize 
the  integral  formula  (6.3). 

Observing  the  recursion  property 


=  (-V-ir)  [  $(Ao;  f,?)  -  ?  ( fzd 


(6.4; 
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we  obtain 

/  i I  ?}  [(&*£  + XQ)  f  (%)[  drt  (p) 

\5Z 

-  -  £  f  <#  *K)  ■$(*..  K,  t?)  +**)  f<ri  *°(?> 

*  /  vl.  ?  lfpU^  *  K)  fyij  d, , (?) . 

Integration  by  parts,  i.e.  application  of  Green's  identity  yields  for 
a  function  <{  e  C(lf)(.Sl ) 


52  1 

a  1  * 

In  the  same  way,  by  integration  by  parts,  we  get 


f  fiti] 

>2  >.->» 

“  "3  (  V  °L*(z)  * 

SI 
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Therefore,  we  have 

/  fr*.;  f, 7)  [(&’,+*.)  fty] 

Si  (°-6> 

-  -£  /  r,?)[ ^  (?) 

-  3  /  “5 ^7} 

ft 

provided  that  ^  is  a  four-times  continuously  differentiable  function 
on  SL  .  Thus,  by  combination  of  (6.3)  and  (6.6),  we  have  for  / e  CwfSl) 


Kt) 


21  f  fty  ^  (f  t)  d*  (f)  (6-7) 

n.*o  Sc 


('IT  l  1a°.K;  f.rflA’S^lAYKW^l 


More  generally,  by  successive  integration  by  parts  we  obtain  in  connection 
with  the  definition  of  ^  (*0  (  ■  ■  ,  f,  y)  and  the  formula  (5.13) 

the  integral  formula: 


Let  be  1  f i xed  poi nt  of  the  uni t  sphere  Si  .  Let  ^  be  a 

( 2m  +  2)  -  times  conti nuously  di fferentiable  functi on  on  Si-  .  Then 


I'D 


2T  ~  f  Pr>  da'?)  (6-8) 

nwo  SL 


+ 


(-  -M 

v  V-ttV 


SI 
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1 


Inserting  the  addition  theorem  of  spherical  harmonics  into  the  first  term 
of  the  right  hand  side  we  find 

f(f)  -  Z1  IT  (f)  /  $».;.(?)  do  (?)  (6.9) 

Thi s  formula  gives  a  comparison  between  the  m-th  partial  sum  of  the  ortho¬ 
gonal  expansion  of  f  i nto  spherical  harmonics  and  the  functional  value 
of  ^  taken  at  the  poi nt  g  £  respectively  (with  exp  1 i ci t 

knowledge  of  the  remai nder  term]. 

The  formula  (6.9)  may  be  considered  as  the  theoretical  background  for 
problems  of  interpolation  and  best  approximation  by  spnerical  splines.  In 
order  zo  see  this  we  have  to  modify  our  integral  formula  (6.9)  using  Green's 
function 


(6.10) 

-  /  ff*. . i.4)  ir> . j-,- cfc/o 

SL 


with  respect  to  the  operator  ( •  •  •  (  ^  and  the 

oa  name  ter  ~  e  $Z  . 


Inserting  (5.10)  into  (6.10)  we  obtain  as  bilinear  expansion 


■<1 I 


* 


=  T 


&k+  1 


(6.11 


Hence,  it  is  envious  that 


i*'\x  °  •  ,  f  •'/  ) 

.  (■  >.)••.  (A*.  -f/'VA. 

v*t  ■  L  ’ 
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Sunmarizing  our  results  we  therefore  obtain  the  theorem: 


Let  be  a_  fixed  point  _of  Q  .  Suppose  that  ■£  ^s  _a 
(2m  +  2)  -  times  continuously  di fferentiable  functi on  defi ned  on  • 
Then 


id) 

=  n  21  /  4<?)  (?)  ota/?) 

«co  j  **  *  Si. 


(6.13) 


4- 


.)f  H, ..  (J„,  f,f  fw  W-W 1  do  V  ■ 


G 
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7.  THE  DIFFERENTIAL  EQUATION  (  A*  )  „  f  **  (*"+*•)- (&”+**)  f  - 
Let  S  be  an  element  of  class  p„  of  the  general  form 


■S  (?) 


Then,  observing  the  differential  equation 


1 7  •  i ) 


( A*  +  ^rv)  ^  (?)  *  0 


for  *a  =  O 


nr T 


and  for  all 


T 


€  SZj 


vie  have 


(AmjLSW  3  (Aj*X0)...(AyX »)  Sty  (7-3) 

-  21  ZT  c,: 

*»“o  _/=»-»  4 

-  o 

for  all  €  Q,  . 

On  the  other  hand,  it  can  be  deduced  from  (6.13)  that  any  solution 
£  £  q  1  *>*»-*■  z)  ^  q  )  of  the  homogeneous  differential  equation 

(&y  (&m7  +  s(?)  =  o  (7.4) 


is  representable  in  the  form 


H 


[7.5) 
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with  coefficients  ^nj.  given  by 


/  &(%)  (?)  dur(p)  . 


Si 


(7.6) 


But  this  means  that  d1**  is  the  null  space  of  the  operator 


(A’)m  -  (A’+K)  ■■■ 


(1.7) 


The  integral  formula  (6.13)  will  be  used  now  to  discuss  the  general 
differential  problem 


•  •  •  (  A*  +•  f(?)  =  %(?)  (7.8) 

for  £  €  a  and  Q  €  C  (  SZ )  . 

From  Green's  surface  identity  it  follows  that 


/  [(  A*  +  Aa)  ...  +  $(%)]  S(%)  do  (7.o) 

-  I  +  •  •  •  (a; +Am)  $(?)]  ftt)  d* 

Si 

=  o 


for  all  elements  S  belonging  to  It  is  clear  that  any  function 

S  £  can  be  added  to  f  without  changing  the  differential  equation 

(7.8).  However,  if  we  require  that  £  is  orthogonal  to  the  null  space  of 
(  a*  +  A0)  ...  (  A'  >Am)  ,  then  the  differential  eauation  is  uniquely 

solvable: 


Jk _ 


A 
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Let  ^  be  _a  continuous  functi on  on Q,  orthogonal  to 


SI 


f  f(?]  o'* 


=  O 


(7.10) 


for  ki  =  Of  ...  f  m.  j  j  -  st  t  ...  t  2.  n  +  a 
Then  the  function  given  by 


'-f  (f)  ~  (~  J  -L; f'£)  doty 

Si 


(7.11) 


represents  the  only  (2m  +  2)  -  times  conti nuously  differentiable  solution 
of  tne  di fferential  equation  (7.3)  which  is  orthogonal  to  ■ 

For  integers  mn.  ^  A  the  bilinear  expansion 


#  ...  ,  •  ffy)  =  (‘fir)  2Z 


Ik*  A 


S(f 


7> 


is  absolutely  and  uniformly  convergent  on  the  unit  sphere  <rz  .  Inter¬ 
changing  sum  and  integral  we  therefore  obtain 


fV?)  - 


H 


Zk  +  4 


k*m+'!  (\~K)  SL 


f  ^  ( f ?)  d* 


as  unique  solution  of  the  differential  equation  (7.8). 
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8.  THE  HILBERT  SPACE  ( SZ) 

In  the  class  C  C  Y<^2)  of  all  (2m  +  2)  -  times  continuously 
differentiable  functions  on  the  sphere  SZ  we  introduce  the  inner  product 

(4,  3-L  <3'1) 


=  £  £(/WJ  v  '//WJ'. w 

The  space  C  (2>~l  *  2  }  ( S2,  )  equipped  with  the  inner  product  f  -t  ■) 
is  a  non-complete  inner  product  space.  For  the  sake  of  simplicity  we  use 
the  abbreviation 

^  ?)  “  ^  *5  $,?)  .  (8.2) 

The  kernel  (^  ,^)  admits  the  series  expansion 


oo  in*-- 1 


(  t, 

-  z:  2 

"r«  v 

^  %  (f) 

Vw- 

(8.3) 

where  the  functions  VM  ; 

n'<f 

are  defined  as 

fol 1 ows 

/  — 
v  2n+-t 

'  ^ 

n.  =*  0;  .... 

/  =  *1  ,  . . .  (  2  n  + 

ll 

>-r 

1 

f -  ^ir) 

rn  +  * 

c 

00 

/  “  * 

/  -  -*  ,  ■  ,  in  *- 1 

• 

s 
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( f  ,  %)  has  the  representation  of  an  isotropic  covariance  function 
on  the  unit  sphere  ^2  (  with  respect  to  Legendre  polynomials).  It  is 
clear  that  if , %)  is  symmetric  in  the  arguments  if  and  p  , 

(ft?)  =  (y ,  fj  •  the  Hecke  formula  we  have 


V-TT 

£n+'t 


(8.6) 


o 


n.  &  k. 


Y-ir 
Zn  +  * 


y\  =■  k-  . 


Consequently,  an  easy  calculation  shows  that 


(  f,  (f,  ■))„  is. 7> 

yyy  -t-'l 

-  ZT  ZT  ■£*,.(%)  f  $(?)  ,  (<7  )  do  ( p 

n.*o  J.*  ‘<f  '  si  '*  C 

+  ( -fc)2nt4i  f  fa; +\)  ■  .  ■  ?#v^)  ■  <n*h). 

si 

Cut  this  means  that 


it)  -  <  4,  (f,  -))^ 


for  every  function 


4  e. 


and  every  point 


(  f,  y  )  is  a  "rep- oducina  kernel"  in  C  <irr<  r‘*'>  )  . 

However,  the  kernel  k"m  if  itself  does  not  belong  to  C  ■;  O. )  -. , 
a  function  of  ^  for  fixed  f  (or  as  a  function  of  If  for  fixed  y  ). 

And  that  is  the  reason  why  we  have  to  complete  our  soace  c? 1,21  (SZ.)  witn 
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respect  to  the  topology  induced  by  (■  ,  •) 
(  Zm  +■  X) 


Let 


all  functions  ^  satisfying 


(Si)  denote  the  space  of 


t-n+4  2 

-  n  niu.r^j  < 

Jm  + 


(8.9) 


Then 


(8.10) 


j_s  the  reproducing  kernel  for  the  Hi  1  bert  space  &  *  ^SZ  )  , 


i  .e. : 


( i )  For  each  fixed  Jf  6  Q.  ,  (fi?)  considered  as  function 

of  £  ,  isin 

( i  i )  for  every  function  $  &  JC  (*~'  **USl )  and  every  point 

$  £  S2  the  reproducing  property 


■fff)  -  (4,  *„(%,■)) 


(8.11) 


is  valid. 


fhe  system  j  Y  l 

- i  «,j,  / 


n.  =■  O,  si  '  . .  . 
J  ‘ 


is  a  Hilbert  basis  in  7f  <im*ZYs c) . 


MOTE:  For  the  case  m  =  0  see  Freeden  (1978/1980)  and  Krarup  (1979). 


In  the  sequel  we  shall  analyze  in  detail  that  a  proper  abstract 
setting  fo*-  the  problems  of  interpolation  by  spherical  splines  will  be 
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provided  by  the  Hilbert  space  SC  dm**-)  (  q  )  under  consideration. 

For  that  purpose,  we  have  to  prove  first  that  SC  C -Sd  ) 

a  Hi ibert  function  space  of  conti nuous  functions  on  the  unit  sphere. 

Let  be  an  element  of  SC  (  SS )  .  Then,  by 

Schwarz's  inequality,  we  get  with  a  suitable  positive  constant 
Cm  ($ ,  f)  dependent  on  m 


/  (Vi  (8.1 

uniformly  with  respect  to  all  points  'f  £  &  .  Furthermore,  there  exists 
a  sequence  (4^1  of  elements  in  C'  ( Si)  with 

^  (in  the  sense  of  the  norm  ( •  ,  •  )rrx)  . 

Therefore,  in  connection  with  (8.12),  we  obtain 


/&(?>  -  I  *  /f  (8.i3) 


uniformly  with  respect  to  all  if  £  .  But  this  implies  that 

+  a  (  'f )  - •-  (  $ )  uni  fonnly  wi  th  respect  to  'f  £  S2  . 

Consequently,  ^  is  a  continuous  function  on  the  sphere  SC  . 

The  Hilbert  space  SS  (  ***  Z>  (  S2.  )  is  naturally  equipped  with 
the  (Sobolev- like)  semi -inner  product  <  •  ,  •  >  defined  by 


SI 


(8.14) 


correspond! ng  to  the  semi -norm 


l£l 


/  -/ 

i  V-T 


)  yfi 

Si 


(3.151 
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where  (A* )  ^  m  (  A"+  A„)  . . .  ;s  to  be  interpreted  in  the 

distributional  sense. 

p 

The  null  space  of  the  semi-norm  II  =  /<•  ■>  is  known  to  be 
simply  the  linear  space  ^  of  dimension  m  =  +  -i  )  A . 

It  should  be  noted  that  the  norm  may  be  physically  interpreted  (at 
least  for  m  -  o  )  as  the  bending  energy  of  a  (thin)  membrane 
spanned  wholly  over  the  unit  sphere  £2  ,  <f  denotes  the  deflection 

normal  to  the  rest  position  of  course  to  be  spherical.  This  model  is 
suggested  by  the  classical  interpretation  of  the  integral 

b 

/  /  d x 

a 


as  the  potential  energy  of  a  statically  deflected  thin  beam  which  indeed 
is  proportional  to  the  integral  taken  over  the  square  of  the  (linearized) 
curvature  of  the  elastica  of  the  beam  (cf.  Freeden  (1981b)). 

We  summarize  our  results  as  follows:  The  (semi-normed)  space 
w  **’(  Q)  defined  by  (8.9)  and  (8 . 15)  a_  semi -Hilbert  space  of 
continuous  functions  on  the  sphere  S2, 
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9.  THE  HILBERT  SPACE  ( Q) 

Let  ? *  i  ■  ■  -  ,  7m  be  a  fundamental  system  of  order  m  on  Q,  .  Then 
there  exists  a  unique  6  £  u'p  satisfying 

£(?*)  -  y*  ,  k  •  (9.D 


for  any  prescribed  (real)  scalars  y^  ,  ■  ,  yM  .  For,  as  function  of 

(p  £  is  representable  in  the  form 

/ 


S(f) 


n  r 


n • o 


f 

f 


?)■ 


(9.2) 


Substitution  with  £  =  ^  Ay ^  gives  M  linear 

equations  in  these  M  coefficients 


c~,\  ( “7  )  ~  V 

(9.3) 

/“  ^  0*)  =  y„  • 


The  linear  equations  are  uniquely  solvable,  since  their  matrix  (cf.  (4. 
is  assumed  to  be  of  full  rank  for  a  fundamental  system  of  order  m. 

r* 

Hence,  given  a  fundamental  system  ^  .  .  ■  ,  of  order  m  on  p 

then  there  are  determined  uniquely  M  linearly  independent  elements  of  ^ 
3,,  ,  .  .  ,  such  that 
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<■*1  ln+4 

71  H 


*a  *-o  J*  -i 


^Vfe  )  ~  dik 


(9.4) 


2«  (?J  ■  *  •  3*  (h)-  0  ,  ■  •  •  /  3  ^)  -  ° 

Si(rj  --  °  ,  3z  (?*.)  =  ^  ,  •  •  •  ,  ^  -  0 


3M  (?J  '  0  -  3,/?*)  “  0  /  •  •  *  ,  3m 

For  any  we  have 

£(?)  =  £  -S<>  )  3fc  C?)  (9-5) 

k.  ~  'I 

as  the  (uniquely  determined)  ^  -  interpolant  of  S  on  the  point  system 

C  4  I  '  '  '  t  *  M 

For  any  choice  of  (real)  values  y ^  t  .  .  yM  ,  the  element 


S(l'  -  p  y„  3hff) 

is  the  unique  solution  of  the  interpolation  problem  £  fy^)  -  yb  t 

k  -  ■  • .  ,  M  /  ^  sir 
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The  polynomials  ,  ...  ,  3 M  of  the  class  ^  are  called  the 
fundamental  polynomials  for  pointwise  interpolation.  In  order  to  evaluate 
the  functions  given  by 

ivt 

=  T_  71  SM/  (fi 

n  =  o  j  * i  v 

corresponding  to  a  prescribed  fundamental  system  ,  • • •  ,  of  order 
m  on  SI  we  have  to  solve  the  linear  system 


*n,  Im-t' i 


(r.) 


\ 


0  \ 

A  i 


0 
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i  .e.  the  matrix 


is  equal  to 


(«”)'*  -  K^Y)'1 


(( 


transposed  matrix). 


For  every 


4  * 


X (SI 


21  4??*)  5*  ( ?) 

hr  * 


is  called  the  (general  i zed)  Lagrange  form.  The  mapping 


Jc  *  * 


is  a  linear  projector  of  X  *^(3. )  with  range  and  kernel 


fin i*2) 


(SI)  -  {feX'^lffv-O^.-f,. 


Eacn  function  e  £(**•+*)  j  can  represented  uniquely 

in  the  form 


<  fo  (?)  , 


(9.10) 


where  (f0  is  of  class  X(*Z***  (  si)  and  £  and  £  ^  ar 

orthogonal  in  the  sense  of  the  inner  product  <  •  ,  •  >ryv : 


A,  4  >, 


(9.11) 


(&)*“* /•  ht* VA')  ■  ■  ■  (*",**")  4  Vl  ^ 

Si 


O  . 


Equipped  with  the  semi -norm 

11 


^  ■  ,  *  >, 


the  linear  space  (Si)  is  a  Hilbert  space  of  continuous  functions  on 

the  unit  sphere  J2  In  view  of  the  reproducing  property  of  the 

Kernel  *  (f,n)  in  7C  +  '  ( SZ  )  it  is  easy  to  see  that 
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C  ,  f  "VY,  ■  ,  ?,  ->m  (5-12) 

-  <  f. ,  f  _ f,  ?j  3„  >m 

n*  a 

-  <  fol  T.  2^11) 

'i 

+  <  £  ,  Z  £  3*1'?;  ? "YV--. -I-. • 3, >m 

n- *  42  =  * 

-  /o  f  ?; 


for  every  function  ^  0  and  every  point 

•f  €  r? . 

Moreover,  the  kernel 


*C  (I,?) 


iw(K,  ?, ?) 


(9.13) 


M 

~  Jz,  "f  ’  (A°, 


-  H  *«'?>  y«YJ. . 


A.=  >» 


A*f  M 


+  Jz  ^  ^  '  *«v  AtC*' 
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belongs  to  the  space  jg  <**«+*■)  ^  . 

Therefore,  our  considerations  can  be  summarized  as  follows: 

The  linear  space  <^f0  **  ( •&)  defined  in  (9.9)  is  equipped 

wi  th  the  norm  f  -  I**,  -**  i/*-  ■  r  a  Hi  Ibert  space.  <**'+*)  fJ2) 
possesses  the  reproducing  kernel  (9.13). 

a 
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10.  OPTIMAL  INTERPOLATION  PROBLEM 


The  basic  problem  of  spherical  interpolation  shall  be  repeated  now 
explicitly  as  follows: 

-  given  an  admissible  system  ^  ,  •  ■  •  ,  ? „  of  points  of  the  unit 

sphere  a  and  a  set  ,  ■  ■  -,y„  of  real  scalars, 

-  construct  a  function  5  .•  Q  — *■  belonging  to  the  linear 

variety 

vn  =  {  f  €  f(7k)  =  yM ,  k- W 

of  all  interpolants  in  ^  )  to  the  data. 

To  assure  uniqueness  of  this  problem,  additional  information  is  clearly 
necessary.  In  problems  of  determining  geodetic  or  geophysical  quantities 
some  restrictions  on  smoothness  and  polynomial  precision  are  required.  This 
can  be  achieved  by  restricting  the  set  of  interpolants  using  an  (energy) 
semi-norm.  An  interpolant  minimizing  such  a  semi-norm  can  be  regarded  as  the 
"smoothest"  solution.  In  this  connection  the  semi -norm  I  •  I  represents 
a  natural  setting  as  explained  in  our  introduction. 

From  the  definition  of  |  •  I  it  is  clear  that 

^  ~  (l^)  /( l(ty')o)  .  3n  (y)jt ota 


o 


For  n  =  '1 


M 


•''ence,  the  minimization  problem  (  N  >  M) 


I  s| 


V'  = 

/v 


<6  v„ 

it"*  1  ,  I 


i  e'  i<?t\ 


i  \  0 . 3 ' 
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is  equivalent  to  the  problem 

l  |s.lm  -  l£9m  j 

if 

j  ^  A/})’  (10'4) 

Remark:  The  case  N  =  M  reduces  to  strict  interpolation  by  spherical 
harmonics . 

From  Chapter  9  we  know  that  /  *z)  1S  the  reproducing 

kernel  of  X  e*"'+t>  ( SI )  .  Thus,  for  f  ,  %  €  SI, 

<  «°n-,  •) ,  *C  (■,?)>„  (!0'5 

-  C (Sffflq+V ■  -  ■  WJAJ *C <■  <’ ?! d “ 

>52 


But  this  means  that  the  matrix 

j 

(10.6) 

is  symmetric  and  positive  definite,  viz.  as  Oram  matrix  of  the  (linearly 
independent)  functions  f ,  p)  ,  •  •  •/  * "m  f  f 1  • 

these  results  now  can  be  used  to  prove  the  main  theorem: 

Suppose  that  ( •? , ,  >; )  ,  .  .  .  are  prescribea  data  points 
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corresponding  to  an  admissible  system  ^  ,  .  .  .  ;  %  ^  of  order  m  on 

the  unit  sphere  .  Then ,  the  functi on  s  €  gi ven  by 

S(? )  w  (?)  +  s°  (?)  (10. 

with 


-V?)  - 


and 

5c  f  7) 

-  (10.9) 


M 

1 z 

n  =  -i 


y«  (?) 


(10.8) 


A*/ 


a.  K 


(  ?*,?) 


is  the  only  solution  of  the  interpolation  problem 


c  n  £ 

4*  v* 


/ 


i .  e . : 


( 

’  _L 

<^7T  / 

.  .  (&"**„)  s(r)\Ldj, 

-  in  4 

/  * 

f^r/ 

1 1 

SI 

where  the 

coefficients  Q„..  ,  .  . 

.  ,  Qm  satisfy 
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the  1 i nea r  equations 


A/ 


-  £  -  £  yn 


(10.10) 


O  *  M  ^  /  ••  •/*/)■ 


Proof.  (Uniqueness)  Suppose  that  s  0  is  given  in  the  form  (10.9). 

Then  is  uniquely  determined,  since  the  linear  equations  (10.10) 

admit  one  and  only  one  solution  >  ^n)  because  of  the  positive 

definiteness  of  (10.6). 

(Minimum  property)  In  order  to  prove  the  minimum  property  of 
in  V%  with  respect  to  I  I  ^  we  consider  the  difference 
of  and  any  function  e  '• 


<•  • 


(10.11) 


Of  course,  we  have  <^0  ^  )  -  O  for  k  =  M  +  -i ,  .  .  .  (  /V  . 

Hence,  by  virtue  of  the  reproducing  property  of  K  in  *l>(Sc) 

it  follows  that 


^  A  J  do 


Si 


(10.12) 


A' 

'  ^ J  ■  ••  h^-)  <4  A.)  <C  at. 


o 
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Therefore  it  is  obvious  that 

^ I  ^ ^7  •  ■  •  ( ^7  foCp)}2,  do3  (p)  { 1 0 . 13) 

=  (£)*"mi/I(a;+a.)  ...  (a;  +  >„)(s0(?)~  d0(?))izda(z) 

si. 

=  (&)****/  |  (&*+\)  •••  s0(?)lz  doty 

SI  7 

*  I  temr+**)  •  •  datylz  doty  . 

SI 

This  shows  that 


U)*™  1 1 (*;+*.) 


(  A.’+*„)  S,tt)lL  duty 


-  (£.)S~'1  J  l(A’  +  \.)  ..  •  {.(yf  duly 

sx> 


for  all  e  1/^ 

T-e.  fa  =  sQ. 


where  equality  holds  if  and  only  if  atc  =  o 

a 


For  later  use  we  write  the  matrix 


U 


►n 


H  +  ^ 


1  X  , 


'  / 
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in  decomposed  form.  Observing  the  explicit  representation  of 
(see  (9.13))  we  obtain 


(  ^0|-|  ;  7*t>i  iVh+a)  *  *  *  ^ 


g"  =  ' 

M 


(10.14) 


?  ^0/  -v  ■■■  'f1 


,  \ 


5A„) 


5  (^o,  •••/^;  /V/fcJ-  ••  /  ?*,  ini  J  \^ri  ^7n^)  ■  ■  ■ 


3 „(&*,)  ...3, 


t ,■&*,/&;  •  •  •  f ‘‘’Vv W^;v 


i % 3M(fa'i  I  \f  7*,, 


51 


Vo! 


3.  (?N)  ■  •• 


3<3*j  •  -3/^ 


\  rVA0)  -,  A  •*;&,£,)  •  •  •  ^o,..,A^,7*j 


VO-tyrJ 


According  to  the  definition  of  the  fundamental  functions  3,,  ,  •  •  •  , 
of  the  class  we  obtain 


••• 


(10.15) 


A  1(7^ 


"t  .£««♦* 


2T  2-  c„  5  •  \ 

~T0  "7 


*»  i*t+i 


Z  z  <4  $..*  <>j 


►t»c  /  «■  >» 


J 


J 
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/ 
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^  ^ 


$o,Ul«) 


c  ( </  'j  \  “  i 


^  $>0Jih)  ••  *  5«****^ 


3 
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11.  SPLINE  FUNCTIONS 


The  consideration  given  in  Chapter  10  now  leads  us  to  the  definition 
of  a  sonerical  spl i ne  with  respect  to  an  admissible  system  of  orc;i"  m  on 

a  : 

Let  be  an  admissible  system  of  order  m  on  the  unit  sphere 

SZ  .  Then  any  function  s  e  #  «*•*♦•*»  ^ of  ^he  form 
(10.7) 

M  K' 

5 tv  =  t  +  ^  .  n 

n»-t  *  *  1  il  •  i  I 

with  arbitrarily  given  reals  ,  ...  ,  yN  and  coefficients  ahJ 

is  called  a  spherical  spline  function  in  7C  ami)  ('Q}  relati  ve  to 

I  •  ■  •  I  f/V  ’ 

The  class  of  all  spline  functions  5  in  (S 2.)  relative 

t0  ,  •  •  w  TV  is  denoted  by  (%  ,  ■  ■  ■  ,  ?„)  . 

The  space  7^  (zif.  •  ,  ?v)  is  a  N  -  dimensional  linear  subspace  of 

7t  iim*v  ( v2 )  containing  the  class  . 

From  Chapter  10  it  is  known  that  the  following  interpolation  property 
is  valid  in  ft  <Xn'*s>  (  Si)  : 

Given  if  €  j if  **  1  ( jZ  )  .  then  there  exists  a_  unique 

element  s  e  K*  ( T*,  -  • fa)  satisfying 

=  {(•?*)  my* 


A 
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for  k.  *  'i  ,  .  .  .V  •  Denote  t hi s  uniquely  determined  el ement 

2!  ■',!,)  &  v 

By  virtue  of  the  definition  of  a  spherical  spline  function 
it  is  easy  to  show  that  the  fi rst  i nteqral  relation 

/  x  iht+i  C  . 

'  <hr)  J  l  (*;+*.)  •••  $(f)\  olC3  ( 

uC 


( </-W  )  f  I  (^1  t  *°)  *  *  ’  ^  yy\)  (  ^^)  ~  0^ 

Si 


*  (  f?)  J  ■*  rfa 


holds  for  every 


/  e  3f  ■ 


Moreover,  it  is  easy  to  see  that  the  second  integral  relation 


( 


>-r  ’ 


JI 


(ffy  ~  ^  (11 
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(  r)  f  l  (A? -r  A  0)  ...  f &*  +  ($(>[)  ~  ^  fy)l  zrtio 

viZ 

+  ^  I  (^t  '  '  ‘  +*».)(*(?  ~S4^)\Z~d^ 


holds  for  every  ^  G  (il**  +a') (St)  and  every  S6  /£,,  .  .  ., 

(cf.  Freeden  1981a). 

a 

In  order  to  calculate  an  interpolating  spline  function 
M  __  a>' 

V?J  *  21  ZZq  kl(r.,7)  I11-4' 

r  *  =  -<  fe=#f+'i  “  <- 


we  have  to  solve  the  linear  equations 


A/ 

H 


k° 


(tk'tj.'i 


M 

-  z: 


*<=■* 


% 


<V 


(ll  -S) 


(f- 


whose  coefficient  -  matrix  is  symmetric  and  positive  definite.  Solving 
the  system  (11.5)  is  therefore  a  most  simple  problem:  tne  matrix  can  be 
factorized  by  the  well-known  Cholesky  procedure.  This  can  be  carried 
through  without  any  need  for  pivoting  or  scaling  (cf.  Cnapt.  13). 
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12.  APPROXIMATION  OF  LINEAR  FUNCTIONALS 


Let  us  consider  linear  functionals  J  ■  PC (i  ~x)(Sl  )  —  1R  0f 
the  following  structure 

°  (12-D 

SI 

where  the  function  g  is  assumed  to  be  piecewise  continuous  on  SI  ,  the 
weights  are  real  and  the  points  ,  .  .  .  ^  lie  on 

dTZ  (d:  positive  integer) . 

Though  this  is  not  the  most  general  class  of  functionals  that  might  be  con¬ 
sidered  on  the  Hilbert  space  (£2.)  it  is  adequate  for  most 

purposes  and  applications.  Obviously  included  as  special  cases  are  the 
follow! ng: 

(i)  the  integral  over  Cl  (c f.  Freeoen  (1978)) 

(  (  f )  d<y- 

SI 

or  any  subdomain  ^  of  SI 

I  f(?)  <f  (?'* 

(ii)  the  orthogonal  coefficients  of  a  function  /  on  St  (cf.  freeden 
(1980/1981  a)) 

( 


(12.1a) 


(12.1b) 


(l?)  CL* 
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being  required  in  connection  with  (Fourier)  series  expansions  of  a  function 
into  spherical  harmonics 


(iii)  any  functional  value  ^ f )  ;  j£  £  Q,/ 


<f(1) 


(I2.1o) 


(iv)  any  finite  linear  combination  of  functional  values  of  ■$.  at 
prescribed  points  & '  •  '  fd  on  SZ  (cf.  Freeden  (1981a)) 


r  \  f<v . 


(12 Je  ) 


Our  purpose  is  to  approximate  a  linear  functional  of  the  form  (12.1)  by  a 
linear  combination  L  of  the  form 


Lf  *  -  £  3.(7  )p7f\ 


where  is  a  given  admissible  system  of  order  m  on  ~>c  . 

C 

The  functional  L  is  called  exact  for  the  degree  m,  if  L  S  =•  J  5 
whenever  •£  £  .  The  remainder,  when  L  is  used  to  approximate  J  , 

is  a  linear  functional  R  defined  by  “R  *  J  -  L  .  If  the  approximation 
of  J  by  L  is  exact  for  the  degree  m,  then  "T2  S  =  O  whenever 


Let  {  be  a  function  of  the  space  <J??  (  5L) 

in  the  form  (9.10) 


given 
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{(*1)  -  ^(?)  +  fo  (?)  , 


(12.3) 


where  is  of  class  and  ^ o  is  an  element  of  3froa^1](  si). 


Then  we  fi rst  get 


/? 


*4. 


(12.4) 


provided  the  approximation  J  by  L  is  exact  for  the  degree  m.  Moreover, 
the  reproducing  property  of  ^  rU  (  ~f  ,  ?  )  in  ^  *z')  (SI) 
implies  that 


*/.  - 


(v-ir)  ^  f [(&+K)  ■  ■  -  (A^wJfC  A„) . . .  d<j(f) 


(12.5) 


f.<r j  ■ 


(Observe:  R  is  a  bounded  linear  functional). 


The  notation 
variable  of 


means  the  linear  functional  R  is  applied  to  the  % 
^  0  (  ~i  >?  )  •  funct'ion  K  given  by 


k'(?)  = 


is° 


($ ,?) 


(12.6) 


I 

i 


j 
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is  called  the  spherical  Peano  kernel  of  order  m  for  the  functional  R. 

Combining  (12.4),  (12.5)  and  (12.6)  we  therefore  find  in  connection  with  (12.2/ 


V 


i'-LV 


J  £  (z)l  da 

3  J 


witn 


(12.7) 


I12-8' 

-  -  L? 

■  -i ~  k':  (?„?)  ■ 

Applying  the  Cauchy  -  Schwarz  inequality  to  the  rignt  hand  side  of  (12.7)  we 
obtain 


I*// 


m+'f  rj* 

/ ^  I  (&px°)  •  ■  •  (Amt+*„)  “(zfdci 


ao 


(12.91 
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We  note  that  the  quantity 


<  K,  K  >  =  f  I  (A*+>0)  .  . .  (a;+AJ  Uf^dLo  (12.10) 

SI 

depends  on  the  nodes  jr  .  .  .  /  y N  and  on  the  operator  7?  ,  but  not 
on  the  function  /  &  1  *"***■'(  QS)  . 

Collecting  our  results  we  obtain  the  following  a  priori  estimate: 
let  J  be  a  linear  functional  of  the  form  ( 12..  1)  and  let  L  _be  any 
approximation  of  the  form  (12.2),  exact  for  the  degree  m .  Then ,  for 
each  function  /  €  X<Zn<*^(SL)  > 

|7/  -  L  /  j  /(J;  J?  -Z  I*  L7  l?)  (  £  y)  (12u) 


The  estimate  (12.11)  enables  us  to  calculate  the  best-approximation 
to  7  ,  i  .e.  the  linear  functional  i-  of  the  form  (12.2) 


Lf  = 


N 

T 


\  L.  3-  (?t 


(12.12) 


exact  for  the  degree  m,  for  which  the  quantity  <  ^  ,  «■  >  m  assumes 
its  minimum. 

The  minimum  can  be  obtained  by  solving  the  uniquely  determined 
quadratic  optimization  problem 
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rJ  aJ 


Jf  l  -  HZI  W,  VCfr„.4 

>  C  tz*t*\+- 1  >  "■  a-M*.,  *»H»*  *- 


yyx  l  n. 


(12.13! 


It  is  easy  to  see  that  the  quadratic  optimization  problem  becomes  its 
minimum  if 


a  k'°  (y.  y  ) 

m  **■ 1  (j.  ,  / 


-  K4 i  f,  ?;> 


(12.14! 


for  j  «  (  a/  .  The  linear  equations  (12.14),  however,  have 

a  unique  solution  in  the  coefficients  aMt-,  <  •  •  -,<2  . 


Summarizing  our  results  we  obtain  the  theorem: 

Let  7  be  a  linear  functi onal  on  f  oc  )  _of  the 

form  (12.1) .  Suppose  that  ( zT^  +  ^  (  .  .  .  a^)  is  the  solution  of  the 

1  inear  equations  (12.14).  Then,  for  each  fc‘  X  r~3'  +'*>  (  S2.)  ,  the 

1 i near  functional  L  gi yen  by 


i  -  2Z  ^  [  f(t  )  -  H  3 


represents  the  best  approximation  to_  J7  ( i_n  tne_  Sense  of  the  ~  z- 
estimate  (12.9)). 


The  approximation  l~  to  0  has  as  a  posteriori  estii 
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/  J/  -Lf\ 


-  i±Y” 


(12.15) 


/  /  (l\^  +  A  j  .  ..  (&W£+\^)  '/fjr)/ ^  dC3 


Si 


The  best  approximation  ("F  -  approximation"  in  the  sense  of  Krarup  (1979)) 
opens  a  new  (non-statistical )  perspective  for  geodetic  purposes  of  prediction 
to  supplement  the  gravity  information,  which  can  be  made  at  only  a  relatively 
few  points  by  determining  further  values  of  gravity  at  other  points  or  to 
compute  best  approximations  of  other  quantities  of  gravity.  Furthermore, 
we  have  a  priori  information.  o 

From  a  theoretical,  but  also  from  a  practical  point  of  view  it  is  of 
great  interest  that  there  is  a  close  connection  between  interpolating 
spline  and  best  approximation.  This  can  be  explained  as  follows  (for  a 
detailed  proof  see  Freeden  ( 1981a), Reuter  (1982)): 

Given  a  function  ^  (C^  'j  .  Denote  by  -S^  the 

uniquely  determined  i nteroclating  spline 


I 


a 


W  I  }l, 


( 12.16) 
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to  tne  linear  variety 


/V 


Then  the  best  approximation  L  to  J 

the  property  that  L  ^  =  J  s 

'  £ 


is  also  uniquely  determined  by 
whenever  /  <£■  %  (  Q^)  . 


In  other  words:  there  are  two  equivalent  ways  to  compute  the  best 
approximation  L  to  J 

(i)  solve  the  linear  equations  (12.14)  (using  Cholesky's  factorization) 

(ii)  apply  the  functional  J  to  the  (optimal)  interpolating  spline  -S.  . 

r 


Remark:  Approximation  of  integrals  based  on  the  idea  of  Weyl's  law  of  uni¬ 
form  distribution  can  be  found  in  Hlawka  (1982). 
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1 3 .  NUMERICAL  METHOD 

Using  generalized  splines  as  e.g.  proposed  in  creeden  (1981  a,b)( 1982)  and 
Meissl  (1981),  one  deals  with  functions  lacking  a  local  support.  Hence  the 
normal  equations  are  full  and  the  size  of  the  systems  is  limited  in  inter¬ 
polation  problems. 

According  to  our  considerations  the  data  relative  to  various  inter¬ 
polation  points  now  will  be  exploited  in  two  steps  based  on  a  combination 
of  well-known  procedures  in  interpolation  theory  (cf.  Davis  (1963)), 
viz.  the  Lagrange  interpolation  and  the  Newton  representation  theorem. 

This  is  of  great  computational  advantage  as  regards  the  number  of  data 
and  the  numerical  effort.  The  standard  algorithms  to  be  used  in  out- 
method  (Cholesky's  decomposition,  forward  substitution)  are  indeed  very 
economical  and  numerically  stable. 

The  great  benefit  of  the  method  proposed  here,  however,  is  that 
solution  process  can  be  formulated  in  a  recursive  way  leading  to  tn 
permanence  property  in  spline  interpolation  problems.  This  will  be 
discussed  now. 

As  we  have  shown,  the  spherical  spline  ■;  «  K  (  (  SL) 

that  interpolates  the  data  points  ( ^  )  ,  •  •  •  ,  ( ?*  /  )  can  be 

represented  in  the  form 


■S(^) 


'f-  x*  \  (z} 

X  *  A 


nr  M  +  r, 


K' 


(7n,?) 


(13.1) 


Assumed 

on 


the  points  ;  .  .  ■ ,  form  a  fundamental  system  of  order 
SI  the  functions  3^  ,  .  .  .t  satisfy  the  property 


S’ 


(  K  3 

V  CR' 


(]')  ON 
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for  n,  te  =  4 1  '  M  (  oT  :  Kronecker  symbol) 

If  /o  *  ,  then 


M 

p(^)  =  ^  )  3^  r^) . 

M  SM 


Given  yA  t  .  .  .  t  yM  ,  the  function  defined  by 


p(t)  -  £.  3,r?) 

«  »yf  L 


is  the  uniquely  determined  solution  of  the  interpolation  problem 


p(rn)  -  yn 


n  -  ^ 


M 


in  the  polynomial  space 


The  expression 


M 

p(?)  -  XL  P(7  )  3„  r?) 


’S  called  (general  i  zed)  Lagrange  formula . 
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As  is  well-known  (cf.  e.g.  Davis  (1963)),  the  Lagrange  formula  has 
one  drawback,  namely  the  lack  of  flexibility  when  passing  from  an  inter- 
polant  to  M  data  to  an  interpolant  of  more  than  M  data.  For  instance,  if 
one  desires  to  pass  from  the  space  to  the  space  'PrrL  ^  by  adjunction  of 
2m  +  3  points,  one  has  to  determine  an  entirely  new  system  of  fundamental 
polynomials  that  are  not  related  in  a  simple  fashion  to  the  old  set. 

In  the  Newton  representation,  however,  the  increase  of  data  can  be 
accomplished  simply  by  adding  of  additional  terms  as  described  now. 


Suppose  that  {  j-  Q  =,  „#1/. 

systems  of  order  m  on  SZ  . 


is  a  sequence  of  admissible 


Then  the  (Gram)  matrix 


7~->  •••  \ 


M+CL 

r 

M 


is  non-singular  for  each  integer  Q.  >  * . 


(13.7) 


^  1  J 


Corresponding  to  the  sequence  { 

sequence  /  L _  „  of  evaluation  functionals  defined  bv 

<-  m  =  -v,  a., . .. 


we  introduce  the 


f  “  { (  ?h  +  q,  )  / 


Q  =  *,*, 


Then,  according  to  the  biorthonormal i ty  theorem  (cf.  Davis  (1963),  Theorem 
2.6.1),  there  are  determined  uniquely  two  triangular  systems  of  constants 
h  with  b-  at  0  such  that  if 
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The  result  of  biorthogonal i zation  can  be  formally  expressed  by  means  of 
determinants 


M+i 


'O 


,Yh^) 


L. 


iy6  (i  .  N 


/f  °  /  I /  \ 

(f»t*  I  l) 


{^v  ?hJ 


U' 


,  l) 


Of  great  interest  is  that  for  the  computation  of  the  last  determinant 
only  the  points  ...  *  are  needed  from  our  prescribed 

admissible  system  of  order  m. 

Consider  a  function  ^  of  the  subspace  spanned  by  the  functions 

u~  <  ,  7  )  ,  •  •  •,  ,?).  i-e-  the  set  of  a11  1inear 

combinations  of  the  functions  u.  0  (*  „)  -  i  • 
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The  representation  (13.9  )  is  called  (general ized)  Newton  formula. 

Remark:  For  the  following  considerations  it  will  help  to  describe  the 
bi orthonormal i ty  theorem  of  Newton  type  in  the  language  of  matrices. 

Let  and  %  designate  the  triangular  matrices  taken  from  the 
coefficient  scheme  (13.8): 


and 


(b. 


1  •  -v  • 

J 


,a 


r 


a 

a 


Then  it  is  characteristic  for  the  bi orthonormality  theorem  (cf.  Davis 
(1963),  chapt.  II,  sect.  2.6)  that 


(2 


<& 


n  ♦  a 

M 


t 

r 


i 


( 13.10a) 


(I:  unit  matrix).  This  is  equivalent  to 


<o 


*4  +Q 


/S 

/ 


{ 13.10b: 


Now.  -S'"  is  a  lower  triangular  matrix  with  non-zero  elements  on  its  prin¬ 
cipal  diagonal  and  is  an  upper  triangular  matrix.  Since  the  matrix 

cr  Mh  * a  is  positive  definite  and  symmetric,  the  decomposition (i  3.  ipp;  is 
actually  a  d  /-*•  -  factorization,  where  is  an  upper  triangular 

and  cf  is  a  positive  diagonal  matrix.  Therefore,  there  exists  a  unique 
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upper  triangular  matrix  r  M  +  <a 
such  that 


with  positive  diagonal  elements 

(  t  *  +  a  1  *  -r  a 

^  M  /  A-f  • 


This  (uniquely  defined)  splitting  of  <3"  ^  *  Gi 
Cholesky  factori zation  (Cholesky  decomposition)  . 


is  known  as  the 


D 


Our  purpose  now  is  to  combine  the  advantages  of  the  Lagrange  formula 
with  the  flexibility  of  the  Newton  formula.  As  a  matter  of  fact,  the  solu¬ 
tion  process  of  optimal  spherical  spline  interpolation  can  be  achieved  in  a 
recursive  form. 


As  shown  in  Chapt.  11,  the  problem  of  determining  an  (optimal) 
interpolating  spline  to  the  data 


^  +  '  yn+a.) 


(13.11) 


is  equivalent  to  the  problem  of  finding  the  function  se  f  <V  *  \ 

**  '  (  M  +  Q  / 

of  the  form 


U7) 


H 


Z Z 

*  *  M+* 


(13.12) 


to  the  vector  a  e  =  (a  ^  ,  ••  •  ,  )  satisfying 

the  linear  equations 


m 
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m  +a 

M 

-  H  x,  <v ) 

ft-s-f  <7  7 

ZZ  **  Or„.r;) 

-  $ 

(13.13) 

j.  =  m  +  -i,  ...  ,  At  +  a 

Using  the  notation  (13.7)  the  linear  equations  can  be  rewritten  in  the 
vectorial  form 


M  +  Q.  'N/ 

&  M  a  =  V  ' 


(13.14) 


where  the  vector  y  &  ,  (  yJ~  (  ,  •  •  •  i  Xn  +  a  1S 

given  by 


M 


%  “  $  " 


*ts-1 


(13.15) 


Furthermore,  setting  for  arbitrary  but  fixed  ^  6  J>d, 


*7?)  s  R®  03.16'. 


we  are  able  to  express  the  linear  combination 


£  *■  , 


a. 


f  7*  ,?) 


(13.17 
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as  scalar  product 


a *  £(?)  . 

M  *f  Gl. 

The  (Q,Q)  -  matrix  <3  M  is  positive  definite  and 

symmetric.  Thus,  according  to  Cholesky's  factorization  theorem,  tnere 
exists  a  non-singular,  upper  triangular  matrix  r  M  +  a  (with 
positive  diagonal  elements)  such  that 


<5 


M  +  Q. 
M 


n wa 

M  / 


(13.19) 


(  (  )t ■  transposed  matrix).  But  this  yields  in  connection  with 
(13.14) 


a 


/  _  *  +  Q 

( 


) 


-1 


-  y  . 


( 13.30) 


Thus  our  expression  (13.17)  takes  the  form 


_  m  *  a  <  /  _  /w  +  <2  ’i-  <• 
/M  /  <  V,  / 


(13.21) 
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Using  the  abbreviations 

yefca,  y*  -  (*Zta)-*  y 

and,  for  each  y  e.  SL  t 

&*(*)€  /  *ffy  m(-rZ+a)~i  U(y) 

we  get 

ft(vn.ay< 

-  (  y’)  fr  £*<>)  • 

Using  coordinates 

y*  «  *a  ,  (V)*  -  fy„!.  ,  •  •  ■  ,  y„*,J 

y  i?1  c  ® ,  fa*/]))*  -  (**,„)?).  •  ■  ,  ^“--q ',ij 

this  can  be  expressed  as  follows 


4  » 


^  7*,  7 ) 


<3 


X. 


r? 
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The  sum  of  the  right  hand  side  of  the  last  formula  represents  a 
formula  of  Newton  type  (for  algorithmic  details  cf.  e.g.  Davis  (1963), 
Meinguet  (1979)).  Furthermore,  the  computation  of  y*  can  be 

organized,  indeed,  easily  by  fontard  substi tution . 

Suimarizing  our  results  we  therefore  obtain  the  following  result: 

Let  ^  '  £ H  be  a_  fundamental  system  of  order  m  on  the 

unit  sphere  SZ  .  Then  each  spline  function  -£  £  JL  ( y .  ...  i 
of  the  form  (13.12)  can  be  represented  in  the  form 


Q 


S(> 


(13.24) 


where  the  vector  ym  £  HQ  ;  ( y  */  -  i  t  , 
solution  of  the  1  inear  system 


iAthe 


/  m  -t-  Q  \  t:  „  '  ^ 

(  th  >  y  *  / 


13.25) 


and,  for  each  but  fixed  y  e  Jc. ,  km  ( y)  (z  HZ 

( k*  (y)  ...  £*  (y))  is  the  solution  of  the  1  inear 


system 


’13.26) 


_  H  +  G 
Herein  c 


/  r-  ‘-1  *  <a  >  - 

M  '  / 

Cholesky  -  factors  of  <5-  Q- 


are  the  uniquely  determined 
i  .e. 
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<5 


hi  +  Q. 
M 


(r 


H+Gt 

M 


'f  r  '1+Q 

/  M 


T  r-i+Q  2$  a_n  upper  triangular. 


Tne  first  sum  of  the  right  hand  side  of  (13.24)  is  a  sum  of 
Lagrange  type,  while  the  second  sum  is  of  Newton  type. 


o 


Remark:  The  Cholesky  factorization  is  one  of  the  best  direct  methods 
for  solving  linear  systems  with  a  positive  definite,  symmetric  matrix. 

The  computer  implementation  of  Cholesky's  decomposition  is  simple. 

Economy  of  storage  can  be  achieved  by  working  only  with  a  linear  array 
of  Q(Q  +  l)/2  elements  consisting  initially  of  the  upper  triangle  of 
+  a  which  is  later  overwritten  with  .  inis 

metnod  is  also  economical  as  regards  the  number  of  arithmetic  operations. 

As  is  well-known,  at  most  Q  square  roots  and  approximately  Q3 /  6  operations 
(1  operation  =  1  multiplication  +  1  addition)  are  required,  to  be  compared 
with  abort  a  3  /  3  operations  for  the  well-known  scheme  of  Sanachiewicz . 

For  the  computationa 1  procedures,  well-known  routines  (in  double 
precision,  if  necessary)  are  available  (cf.  e.g.  Wilkinson  -  Reinsch  (1971)1. 
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The  spline  of  the  form  (13.24)  possesses  the  permanence  property,  i.e. 
the  transition  from  the  spline  s<M*Q)  e  (?■,/ ■■  )  interpolating 

the  data  points  ,  ■  ■  ■  ,(  yh+a)  to  the  spline 

s,M+Q-*/')  e  ■■■  ,'?h+q+*)  interpolating  the  data  points 

(7*.y.)  ,  y  )  necessitates  merely  the  addition  of  one  more 

term  to  the  expansion  of  s  ,  all  the  terms  obtained  formerly  re¬ 

maining  unchanged: 


s. 


(  hi  +  Q  4 


Jr.  y-3-'nl  -  £  >CV  T) 


■<n; 


V) 


X 


N-fQ  +  -1 


<"  ■ 


The  price  to  be  paid  is  the  change  of  the  basis  system  for  the  spherical 
harmonics  of  order  &  m  and  the  biorthonormal ization  process. 


The  reasons  for  the  convenience  of  the  permanence  property  are  the 
'ollowi ng: 


by  virtue  of  the  Cholesky  -  decompos’ tion  of  the  matrix 


<3-r+a  - 


the  i-th  row  of 


•¥  a 

H  ♦  Q 


depends  only  on  “he  rows  j.  *  c  of  the  matrix 


(ii)  the  quantities  /*  ;  km  ( ?)  are  computed  by  forward 
substi tution . 


78 


From  an  algorithmic  point  of  view,  the  transition  from  tne  spline 

M  <2 


H  X,  A,  >r<  r  f  x,'  "31, 

/t  ./  *  -*  0  a 


interpolating  the  data  points 


,  x,)  , 


/  ^  <•  <2.  '  Vm  +  Q; 


to  the  spline 


S  (?)  = 


interpolating  the  data  points 


Qt~ 


K  3n  i?)  <-  Z=  y„"  i  *  (?> 


^  <  x, )  , 


H  +  C  fT  /  '  H  +  <2  +~r 


Q.  t-r  ) 


(T:  positive  integer)  necessitates  merely  the  computation  of  the  additional 
term:,  i  n  O 


from 


f-.t  d 


and  the  continuation  of  the  Cholesky  factorization 


to  o 


C2  +  T 


According  to  the  Pi orthonormal i zation  process,  exclusively  forward 
substitution  (cf.  (13.25),  (13.26))  is  required  for  the  calculation  of 
the  sum 


QrT 

y  ' 

,  -  Q 


W  • 


79 


Remark:  Generalizations  to  heterogeneous  data  (in  form  of  (bounded) 
linear  functionals)  will  be  given  in  Reuter  (1982)  and  Freeden/Reuter  (1982). 

Analogous  considerations  establish  the  permanence  property  in  best 
approximation  problems.  c 


Remark : 


In  modified  form  the  interpolation  method  analyzed  here  has  been 
proposed  first  by  Meinguet  (1979)  in  connection  with  the  Beppo  -  Levi 
soace  9 )  of  continuous  linear  functionals  defined  on  the 

class  of  infinitely  differentiable  functions  with  compact  support  in 
9  -  dimensional  Euclidean  space  Ik  q  ,  for  which  all  the  Dartial 

derivatives 

...  ^  J  -  ^  4.  ...  * 


(in  the  di stri buti ona 1  sense)  of  total  order  m  are  square  -  integrable 
in  1R  ^  .  In  this  case  ,  TC  is  naturally  equipped  with 

the  semi -inner  product  corresDondi nq  to  the  norm 


d 


(dV:  volume  element) 


wnere  every  partial  derivative  is  to  be  interpreted  in  the  distributional 


sense . 


14.  COMPUTATIONAL  ASPECTS 


The  results  developed  in  this  paper  give  rise  to  the  following 
algori thm: 

Seep  1:  compute  the  symmetric  matrix  (N  =  M  +  Q) 


Step  2:  print  the  matrix 

/  So,<  (?*) 
j  <?,) 

I 

!  4a  <7*) 

I 

i 

i 

i 

! 

!  S~,-*  4) 

\  * 

\ 

\  ^  ?*  ') 


's».-  \ 
^*.3  (  7n)  I 

.  I 

I 

f 

I 

I 

4*  4/v)  I 


L 
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Step  3:  compute  the  matrix 


\3h(7^)  •  -  •  (7„)j 


1 ^  <t,)  •••  s°,Jr„)  V 


V  ^m.ti****  )  •  *  *  (7h^j 


by  solving  the  linear  system 


So„(7„)  \ 


\  1 


-1 


\  2r»» (Taj 


I  ■  •  •  3,(?n)  \ 


\  3M  (?m*0  •  •  •  $„(?„) 


Step  4:  compute  the  symmetric  matrix. 


H 

(Mh  I 


I  • 


where 


k'  c 

Wx 


7; 


is  qiven  by  ( 10.14) 


M 
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(  t,  7;) 


C£)(^o,  ... ,  Am;  7i/^) 


M 


VZ)(^o,  ■■■,*-*;  Ti,7k^k(r^) 


21  \(ft)  $U)  (X°r  -  ,x~;  ,jr) 


H  M 


Jz  Jz  3*-{V  Vl)V*r~,x«;7kt?H)3('rj) 


Step  5:  Cholesky  -  factorization  of  the  matrix  (  Km(ti,7-)) 


'U  »  M+*I  . .  .  (  /V 
H  +  l,  ...,  N 


(i')‘  r" 

'  N 1  /  M 


Step  6:  determine  the  vector  y  given  by 

=  ^  y  3  f  ^ 

'  W  *  *  /►?  K7  '  (  m  +  * ) 

H  =  'i  ' 


M 


y, 


A/ 


Xv  -  TL  y„  (yN) 


Step  7:  solve  the  linear  system  (forward  substitution) 


«)*> 


y 


r~) 

Step  8:  choose  tne  print  v  £  /  for  which  tne  interpolating 

spline  has  to  be  computed 
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Step  9:  compute  the  vector  3(?)  ,  (3(^))±  =  (3^(y)  (  ...  t  3n(y))/ 
by  solving  the  linear  system 


/  $>,+  (%)  •  -  •  ^ 

=x 

M?) 

•  • 

y  ,im+4  ^?«)  y 

SteplO:  compute  the  vector  K  (  %)  given  by 

{£(?))*  -  (*Z( /  •  •  •  ,  (l„,  ?)) 

Step  11:  solve  the  linear  system  (forward  substitution) 

(-:)*  k'(7)  - 1<^) 


Step  12: 


compute  the  interpolating  spline  s  of  the  form  (13.24)  at  the 
Doi  nt  y  £  -Q 


Step  13:  continue  with  Step  3 


REMARK:  The  matrix 


(7*)  •  •  •  \  ~  ^ 


i 


is  not  used  anywhere  in  our  algorithm. 


15.  A  SIMPLE  EXAMPLE 


We  consider  the  potential  of  the (.buried! mass  point 


U.(£) 


I?  -  yi 


l£  &  ,  l  0/  0/0.9). 


(15.1) 


Our  aim  is  to  approximate  u  by  a  spherical  spline  function 
interpolating  in  given  points  i  •  •  •  t  ?*/  the  urnt  sPr,ere- 

We  chose  62  points  generated  by  regular  polyhebra 


B(0,±T,±1),  6 ( ±1 , 0 , ±t) ,  g(±T,±1,0), 

JT:  y  (  +  P  ,  ±T  ,  0  )  ,  YtO/^S/iT),  yUt,0,±B),  y  (  *  1  ,  ±  1  ,  ±  1 )  , 
(±1,0,0),  (0,±1,0),  (0,0, ±1), 


where  ,  (3  ,  ^ 


are  given  by 
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We  assume  polynomial  precision  up  to  m 


Figure  1  gives  a  graphical  impression 
of  the  accuracy  for  a  meridian  on 
the  northern  hemisphere  passing 
through  the  north  pole. 


In  the  following  we  have  used  two  times  eight  additional  nodes  to  eliminate 
the  derivatives  between  the  interpolating  spline  and  the  function  u  in  a 
local  area  around  the  north  pole. 


Figure  2  illustrates  the  distri¬ 
bution  of  the  points  of  the 
northern  hemisphere  (projected  into 
the  equatorial  plane) . 

"x"  denotes  the  points  of  T  , 
while  "o"  and  V  stands  for  the 
nodes  of  and  respectively. 


±0.2500,  0.9354) 
±0.3536,  0.8660) 


(±0.1721,  ±0.0713,  0.9825) 
(±0.0713,  ±0.1721,  0.9825) 
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For  comparison  Figure  3  shows  the  graph  of  the  interpolating  spline  with 
respect  to  the  70  points  of  ^  w  ^  t  while  Figure  4  is  based  on  all 
78  points. 


Figure  3  Figure  4 

It  turns  out  that  we  get  an  acceptable  approximation  in  the  neighborhood 
of  the  north  pole  avoiding  a  deterioration  of  the  accuracy  in  all  other 
parts  of  the  northern  hemi sphere. 

Unfortunately  it  appears  that  the  Green  function  fw(A0l...,Am; 
with  respect  to  the  operator  cannot  be  expressed  by  an  plementary 

function.  Thus  the  Green  function  has  to  be  replaced  by  an  appropriate 
expression  which  is  convenient  for  computational  purposes. 

In  our  example,  satisfactory  results  can  be  achieved  based  on  finite 
partial  sums  of  the  bilinear  expansion  (6.11). 

Other  suggestions  for  replacing  Green's  (kernel)  functions  by  closed 
forms  convenient  for  computation  will  be  discussed  in  Reuter  (1932)  and 
C>eeden/Reuter  (1982). 
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